function vspec(varargin)
% plot spectrum of current traces --- 
% 29 Sept. 2003 P. Manis
%
% updated 17 Sept. 2005. 
% new parameters for welch.
% twin is a two-element array indicateing the beginning and ending times
% for the spectral calculation.

global CURRENT VOLTAGE DFILE

rate = DFILE.nr_channel(1)*DFILE.rate(1);
fs = 1000000/rate; % in Hertz.
if(nargin == 0 || isempty(varargin{1}))
    hf = findobj('Tag', 'I_Spec');
    if(isempty(hf))
        hf=figure('Tag', 'I_Spec');
    else
        figure(hf);
        clf;
    end;
else
    axes(varargin{1});
end;
t0 = 1;
t1 = size(CURRENT,2);
r = 0.001*DFILE.rate(1)*DFILE.nr_channel(1);
if(nargin >= 2)
    twin = varargin{2}; % get the argument.
    t0 = floor(twin(1)/r);
    t1 = floor(twin(2)/r);
end;
t=[0:rate:rate*(t1-1)]/1000;


QueMessage('vspec - Finding spikes');
[first_spike, first_isi, nr_spikes, spike_train]=find_spikes2(DFILE, VOLTAGE, t0, t1, 0);

vavs = [];
laves = [];
k = 0;
u=find(nr_spikes == 0); % only for traces with NO spikes...
Hs = spectrum.welch('blackman', 1024, 25);
subplot('Position', [0.1, 0.65, 0.8, 0.3]);
for i = u
%        [vsp, vfr] = pwelch(VOLTAGE(i,t0:t1)-mean(VOLTAGE(i,t0:t1)), 1000, 0, 5000, fs);
        v=VOLTAGE(i,t0:t1)-mean(VOLTAGE(i,t0:t1));
        v=winfilt(v, 0, fs/2, rate);
        plot(t, v);
        hold on;
        Hpsd = psd(Hs,v,'Fs',fs);
%        [lfr lsp] = lombpsd2(fs, VOLTAGE(i, t0:t1));
        vsp = Hpsd.Data;
        vfr = Hpsd.Frequencies;
        if(k == 0)
            vavs = vsp;
%            lavs = lsp;
        else
            vavs = vavs + vsp;
 %           lavs = lavs + lsp;
        end;
        k = k + 1;
        % do a lomb's psd also.

end;
vavs = vavs/length(u);
%lavs = lavs/length(u);

subplot('Position', [0.1, 0.1, 0.8, 0.5]);
loglog(vfr(1:end), vavs(1:end), 'b-');
hold on
%loglog(lfr(1:end)*1000, lavs(1:end), 'r');


return;


